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Abstract 

We present an application of Wertheim's Thermodynamic Perturbation 
Theory (TPT1) to a simple coarse grained model made of flexibly bonded 
Lennard-Jones monomers. We use both the Reference Hyper-Netted-Chain 
(RHNC) and Mean Spherical approximation (MSA) integral equation theories 
to describe the properties of the reference fluid. The equation of state, the 
density dependence of the excess chemical potential, and the critical points 
of the liquid-vapor transition are compared with simulation results and good 
agreement is found. The RHNC version is somewhat more accurate, while 
the MSA version has the advantage of being almost analytic. We analyze the 
scaling behavior of the critical point of chain fluids according to TPT1 and 
find it to reproduce the mean field exponents: The critical monomer density 
is predicted to vanish as n -1 / 2 upon increasing the chain length n while the 
critical temperature is predicted to reach an asymptotic finite temperature 
that is attained as n -1 / 2 . The predicted asymptotic finite critical tempera- 
ture obtained from the RHNC and MSA versions of TPT1 is found to be in 
good agreement with the G point of our polymer model as obtained from the 
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temperature dependence of the single chain conformations. 



I. INTRODUCTION 

For a long time, the prediction of the equation of state of polymers from first principles 
was such a difficult task that it could not be solved without the need of very idealized 
models with more or less meaningful empirical parameters. Perhaps the most successful of 
the early attempts was that of Prigogine et a/.i who considered a lattice and introduced an 
ad-hoc empirical parameter known as the "number of degrees of freedom per monomer". 
Later on, Flory et al. extended this theory to the continuum at the cost of introducing 
some additional parameters, leading to what is now known as the FOVE theory. Another 
approach to the problem is based on the polymer+solvent and polymer+vacuum analogy. 
In this way, the well known Flory-Huggins§H (FH) and Sanchez-Lacombel equations of state 
have also been employed to describe the behavior of pure fluids. 

On the other hand, the approach from liquid state theory has taken much longer to yield 
useful results but has now reached a point where a rather satisfactory description of the 
equation of state of idealized polymer models in the continuum is affordable without the 
need of any empirical parameters whatsoever. The most popular approaches are the Polymer 
Reference Interaction Site Model (PRISM) of Curro and Schweizer,! the Generalized Dimer 
Flory theory (GDF) of Honnell and Hall! and the Thermodynamic Perturbation Theory of 
Wertheim (TPT1).0 The latter has been widely used not only because it yields results that 
are of similar or superior quality than other alternatives, but because it is the simplest and 
more tractable of all them and demands a minimum of information. 



Originally, this theory was developed to consider fluids of associating hard spheres.! 




Later on, it was realized that in the limit of infinite associating strength, a polydisperse 
mixture of polymers was recovered.0 Chapman et al. extended the theory to monodisperse 
polymer fluids and rewrote the equations in a very convenient notation.13 Later, by adding 
a mean field perturbative contribution, Jackson et al. showed that the theory was able to 
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describe the behavior of very different sorts of real fluids.0 Meanwhile, it was realized that 
the theory could be used just as well to consider polymers made of attractive monomer S.E3HH 
In this way, it is possible to describe real fluids with improved accuracy and no need for 
a mean field-like perturbation. Since then, the theory has achieved enormous popularity 
and has been applied to describe polymers of tangent Lennard- Jones beadsEM and square 
wells,il'0 as well as to describe real fluids in chemical engineering applications.iHil Fur- 
thermore, modifications of the original theory have been proposed that allow to describe 
realistic polymer models without the need of empirical parameters.i!i! 

A very interesting issue both from the practical and theoretical point of view is the 
behavior of the critical point of polymer fluids as the number of monomers increases. By 
invoking the polymer +solvent and polymer +vacuum analogy, one would expect from the 
FH theory that the pure polymer fluid should reach an asymptotic critical temperature, 
whereas the critical mass density should become vanishingly small. However, in a recent 
paper, Chatterjee and Schweizer0 have pointed out that this analogy cannot be taken for 
granted because the FH scaling predictions are determined by imposing equal chemical 
potentials, whereas the critical point of pure polymer fluids is related to a phase equilibria 
that results from the condition of equal pressure at a given temper at ure.0 

On the other hand, based on rather limited amount of data, several empirical cor- 
relations have been employed to predict the critical properties of substances such as 
polymethylene.il~H These correlations have predicted widely different behavior, ranging 
from infinite critical temperature^ to finite asymptotic critical mass densityBSS More 
soundly based equations such as the FOVE have been recently employed to support the idea 
that the critical mass density could reach an asymptotic constant value.il However, such an 
approach relies on extending the applicability of the FOVE theory to densities below which 
it was not meant to be used J 

Renormalization group calculations, however, show that the FH mean field theory yields 
the correct asymptotic dependence in the long chain length albeit this asymptotic 

behavior is only reached for extremely long chains and the corrections to scaling are enormous 



even for typical chain length considered experimentally.il 

More recently, new experimental techniques have allowed to measure the critical points of 
longer n-alkanes whose critical temperature lays above the point where thermal decomposi- 
tion starts.EiH^I These experiments show that the critical mass density reaches a maximum 
and then starts to decreased Simulation results of both realistic alkanes^ and idealized 
polymer modelsi!'@ give support to this finding. 

Surprisingly very little attention has been devoted to the study of this problem from 
the point of view of the modern theories such as PRISM, GFD or TPT1. In a recent 
paper, the PRISM was employed in an attempt to solve the question in the framework 
of an analytical tractable theory, leading to different predictions depending on the closure 
employed to solve the PRISM equation.0 Other previous studies using TPT1 plus a mean 
field attractive contribution suggested that the critical mass density should vanish in the 
infinite chain length limit .0 However, such a conclusion relied on the assumption that the 
mean field contribution increases linearly with the chain length, a point which at present 
cannot be taken for granted. Nikitin et a/.@'0 have recently presented a theory that may be 
considered as the simplest possible approximation to Wertheim's theory and arrive to the 
same conclusion as reference^!. 

In this paper we will try to reach further understanding of this issue. In the next section, 
we will review the fundamentals of Wertheim's theory by using somewhat different but more 
physically appealing arguments suggested by Zhou and Stell.E3 We will then show by means 
of general arguments how the scaling laws for the critical properties are closely related to 
the virial coefficients (section III). We will then apply TPT1 to a polymer model and briefly 
describe how to implement the theory in section IV and in the Appendixes. Section V 
describes the details of the simulations performed to test the theory. We then devote section 
VI to present our results and close with a brief conclusion. 



II. EQUATION OF STATE OF POLYMERS 
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A. Preliminary definitions 



We consider a monodisperse fluid of polymers made up of n monomers each. Non bonded 
monomers of the same polymer and monomers belonging to different molecules interact 
through some pair potential, uo, while adjacent monomers of the same molecule have an 
additional bonding potential, <£>, responsible for the connectivity of the polymer. This po- 
tential is such that its action vanishes beyond some well defined inter-atomic distance. If 
eventually, two adjacent monomers were found outside this range, they would no longer be 
bonded and the n-mer would be considered to have broken. In practice, this can be avoided 
by making the well depth of the associating potential infinitely large. 

Alternatively, we may consider an associating multicomponent mixture of n different 
monomeric species, A, B, C, etc. Each of these species interacts with members of its 
own species and with members of the remaining species by means of Uq. Furthermore, the 
bonding potential $ is responsible for associating reactions between monomers of type A 
with monomers of type B, monomers of type B with monomers of type C and so on. More 
concisely, the association reaction taking place is of the form: 

A + B + C + ---^ ABC--- (1) 

As the bonding potential is pairwise, the only requirement that is needed to define the 
: 'ABC ■ ■ ■' complex as an n-mer is that each of the adjacent pairs be found within the range 
of the bonding potential. Clearly, for an equimolar composition of such a mixture in the 
limit of complete association, we are lead to a system identical to that described in the 
preceding paragraph. From a physical point of view, this limit is reached for infinite well 
depths. However, in what follows it will prove useful to consider that the well depth has 
some arbitrary finite value. The system is then made of a mixture of free monomers and 
n-mers whose composition depends on the nature of the bonding potential. 

Considering the similarity between the two systems described, we may get an approx- 
imation for the equation of state for the chain molecule by studying the behavior of the 
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associating system. In order to do so, we will obtain an expression for the free energy in 
terms of the degree of association of the mixture. We will then relate the degree of as- 
sociation to the structure of the system and introduce some simplifying assumptions for 
the n-body correlations. Finally we will take the limit of complete association and get an 
equation of state for the chain fluid. 

B. The association reaction 

Let us consider an associating system as that described before, which is initially prepared 
by mixing in equal proportions the pure monomers, so that there are N monomers of each 
species inside a volume V and the resulting number density of each of the species is p. The 
system will eventually reach a state of equilibrium, whereby a fraction of the monomers of 
each species has associated to form n-mers. Let this fraction be a. Then, the remaining 
concentration of non bonded monomers of each species is given by p(l — a), while, according 
to the stoichiometry of the reaction, the concentration of n-mers will be pa (note that in 
the limit of complete association p will actually designate the polymer number density). 
Schematically, the process can be described as follows: 

A + B + C ■ ■ ■ # ABC ■ ■ ■ 
initially p p p (2) 

at equilibrium p{\ — a) p(l — a) p(l — a) pa 

Obviously, the number of n-mers formed is V pa = Na, while the number of remaining 
free monomers of each species is N(l — a). The total Helmholtz free energy of the system, 
A = G — pV is therefore given by: 

n 

A ( a ) = Ar ( 1 ~ a )^i( a ) + Nap n _ mcr (a) - p(a)V (3) 
i=i 

In what follows, we associate a number from 1 to n to each of the n species. Therefore, 
Pi stands for the chemical potential of species % in the mixture of composition given by a. 
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It is important to realize that the composition of the mixture as given by the association 
degree is appropriate to a given model potential. In the above expression for the free energy, 
we have considered a system of monomers that interact through some reference potential, 
Mo, while a bonding potential is responsible for the connectivity of the n-mer. Let us now 
consider a reference fluid made of monomers that interact through the reference potential but 
that have no association potential whatsoever. In this case, the free energy is best expressed 
in terms of the chemical potential of the monomers alone. For the sake of simplicity, we will 
denote this free energy as A(a = 0), albeit from a strictly geometrical point of view, n-mers 
could form eventually even in the absence of a bonding potential: 

n 

A(a = 0) = J2 N^i(a = 0) - p(a = 0)V (4) 
i=i 

The free energy of the associating system measured relative to that of the reference fluid 
is then given by: 

AA n n V 

-jr = J2l^i( a ) - Vi(a = 0)] + a[fi n ^ mci (a) - ^^(a)] - \p(a) - p(a = 0)]— (5) 

i=l i=l 

This equation can be further simplified by invoking the condition of chemical equilibrium of 
the associating mixture, which reads: 

n 

- £*n-mer(a) = (6) 

i=l 

Substitution of the above equation into Eq. [| yields: 

AA n V 

— = X>*(«) - = 0)] - [p{a) - p{a = 0)]- (7) 

Let us now express the chemical potentials of each of the components in terms of an 
ideal and an excess contribution: 

tn(a) = /4 d (a) + tf(a) =rf + kT\n Pl (a) + ^ x («) (8) 

Introducing the above expression for the chemical potential into Eq. [7], we get: 

^ = ^InA ff,., + f>r(«) - = 0)] - [ P (a)-p(a = 0)£ (9) 



If we now recall that Pi{a) = p(l — a) and pi(a = 0) = p, we finally obtain: 

A A n V 

— = nkT\n(l -a) + EK X («) - ^> = 0)] - ]p(a) - p(a = 0)]- (10) 

i=l 

This is an exact equation for the difference in free energy between the associating system 
(with aN n-mers) and the reference non-associating fluid of free monomers. As it stands, 
it may seem rather useless, because it is a function of the unknown quantities: a, p^{a) 
and p(a). However, we shall see that in the low density limit, the above equation becomes 
a function of the degree of association only and that this, in turn, may be obtained from 
knowledge of n-body correlations in the fluid. In a further approximation, the n-body 
correlations of the fluid will be expressed in terms of n-body correlations of a reference 
fluid of non-bonded monomers. Finally, by invoking a superposition approximation, we will 
express the n-body correlations in terms of two body correlations and obtain an expression 
for the pressure that depends solely on known quantities of a reference fluid of spherical 
particles. 



C. The association reaction in the limit of low density 

In the limit of low density, the equation of state of the associating system will depend 
only on the total number of particles in the fluid, N(a): 

w = N{a) 

N(a) is simply obtained by summing the number of particles of each species: 

n 

N(a) =^N(l-a) + iVa = niV(l-a) + iVa (12) 
i=i 

On the other hand, the number of particles of the completely un-associated system is 
simply N(a = 0) = nN. Therefore, the difference in pressure between the non-bonded 
system and the system with degree of association a is: 

Ap(a)^ = -kTa(n - 1) (13) 



By using this expression and considering that, by definition, the excess chemical potentials 
vanish in the limit of low density, Eq. |1(] becomes: 

AA 

— = nln(l-a) + a(n-l) (14) 

This is an exact equation for the difference in free energy in the limit of low density. Applying 
the standard thermodynamic relationship connecting the pressure with the free energy yields 
the exact expression for the difference in pressure in the limit of low density: 

Ap ail — n) — Ida . 

T^ = P i — 7T 15 

kTp 1 — a op 

In what follows, we will consider this expression to be valid in all the density range. The 
search of an equation of state for the associating system will be thus accomplished if we find 
an expression for the association degree in terms of known properties. Once this relation 
has been found, we will obtain an equation of state for the system of n-mers by taking the 
limit of complete association, i.e., a — 1. 

D. Relation between the degree of association and the structure of the fluid 

1. Expression for the degree of association in terms of the excess chemical potential of the 

components 

First, consider the equilibrium constant of the reaction, defined as the ratio of the con- 
centration of the products to that of the reactants: 

K= , ap , (16) 
(1 - a) n p n 

The connection of the equilibrium constant to the thermodynamics of the process may 
be obtained by expressing the chemical potential of each of the components as in Eq. |8] and 
substituting into the condition of chemical equilibrium, Eq. ^|. After some simple algebraic 
manipulations, we are lead to the following expression for the equilibrium constant: 

n n 

kT\nK + //°_ mer - 5>° + // n x _ me » -5>f(a) = (17) 

i=l i=l 



This expression may be further simplified by considering that, in the limit of low densities, 
the excess chemical potentials vanish. As a consequence of this, the low density equilibrium 
constant, K , is given as follows: 

n 

kTlnK = J2»i-Vn- m er (18) 

i=l 

Substitution of this expression into Eq. [17|, leads finally to a simple equation for the equilib- 
rium constant in terms of the excess chemical potential of the components of the mixture: 

K n 

kT ^ JT = E/4» - ^- mc » (19) 

K o i= i 

2. Expression for the structure of the fluid in terms of the excess chemical potential of the 

components 

In order to relate the structure of the fluid to the excess chemical potential of the com- 
ponents of the mixture, (i.e., A, B, C, etc. monomers and n-mers), let us consider the 
thermodynamic cycle of figure 1. 

In the first step of the cycle, an isolated n-mer is dissolved into a fluid mixture with 
association degree a. Initially, the total Gibbs free energy of the system is the sum of the 
free energy of the isolated n-mer, G n _ mer and the free energy of the mixture, G m i x (a). After 
dissolving the n-mer, the resulting free energy is that of the original mixture with an extra 
n-mer, G mix+n _ mer . The change in G is therefore: 

AG\ Cmix+n — mer ^mix G n — mer (20) 

In the thermodynamic limit, the difference G m ix+n-mer — G mix becomes equal to the chemical 
potential of the n-mer in the mixture, while G n _ mer may be considered to be the chemical 
potential of the isolated n-mer (i.e., the free energy difference between a system with a single 
n-mer and an empty system). It is thus seen that: 

AGi = /i e n x - mer («) (21) 
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In a second step, n uncorrelated monomers of type A, B, C, etc. dissolved in the mixture 
are brought together in such a way that the resulting ABC ■ ■ ■ complex forms one of the 
many possible conformers of the n-mer, say, one such that the vector joining B to A is 1*12 , 
that joining G to B is r23, etc. This event will occur according to a probability density 
given by the n-body correlation function of the mixture, p n g(r 12 , r 2 3, . . . , r n _ l n ). The n-mer 
density is given as an integral of this function over all the conformations compatible with 
the monomer: 

Pn-mcr = p™ / ■ ■ • / 9(^12, r 23 , r„-_i >n )d 3 ri 2 d 3 r23 • • • d 3 r„_i jn (22) 

J V J V 

where v is the volume within the range of the bonding potential, i.e., any two adjacent 
monomers whose distance vector is not within this volume are not considered to be bonded. 

The process of forming the n-mer from a set of n uncorrelated monomers may be con- 
sidered as a chemical reaction of the form 

n uncorrelated monomers n correlated monomers 

Friedmanil has shown that such an equation is characterized by an equilibrium constant of 
the form K eq = p n _ mer /p n . Substitution of Eq. ^ into the expression for the equilibrium 
constant, yields K eq = p n A/p n , where, 

A = / • ■ ■ / g(r 12 , r 23 , r„_i ;n )d 3 ri 2 d 3 r23 • • • d 3 r n _i in (23) 

Now, the free energy of the whole process may be written down as: 

AG = AG 2 + kT In K cq (24) 

However, when the system reaches equilibrium, AG = 0, so that we finally obtain: 

AG 2 = -£;Tln A (25) 

Similar arguments as those put through for the first and second steps of the cycle, lead 
to the conclusion that 
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AG 3 = £^ X (a) (26) 

i=i 

AG 4 = -/cTln A (27) 

where A is the integral of Eq. evaluated at zero density. 

Substitution of Eq. |21], |26], |27| into the net energy balance of the cycle, AGi + AG 4 = 
AG 3 + AG 2 , leads to the desired equation relating the structure of the fluid with the excess 
chemical potential of the components of the mixture: 

fcTln ^ = E^ X («)-K- mer («) (28) 

A similar equation has been derived in a rather more formal way recently^! for the particular 
case of infinitely short ranged association potentials. The derivation we have employed is 
based on the thermodynamic cycle presented by Zhou and StellS, which allows to extend 
their previous result to finite range potentials. 

3. Expression for the degree of association in terms of the structure of the fluid 

Substitution of Eq. [19] into Eq. |28] shows that the equilibrium constant is related to the 
n-body correlation function of the associating mixture through the following relation: 

Finally, using the expression for K in terms of the degree of association, Eq. [Tj], we obtain 
the desired equation relating the degree of association with the structure of the system: 

t, n r = 30 

(1 - a) n p n - % A 

E. The equation of state 

Previously, we obtained an approximate equation for the pressure in terms of the as- 
sociation degree of the mixture, Eq. [1^. The density derivative of the association degree, 
required in such an expression is obtained from Eq. ^: 

12 



da _ 1 - a {n - I) a + pa^ 
dp p 1 + a(n — 1) 



(31) 



Substitution of this result into Eq. [15] yields an expression for the change in pressure due 
to the formation of the n-mers from a fluid of non-associated monomers. This expression 
depends solely on the degree of association and the n-body correlation function of the asso- 
ciating system: 

_JL = _ a(n _ 1+p _) (32) 

By taking the limit of infinite association, which physically corresponds to infinitely 
increasing the well depth of the bonding potential, we would arrive at an equation for the 
pressure of the n-mer fluid relative to that of the monomer reference fluid. In order to do so, 
however, we would require the n-body correlation function of the associating system, which 
enters through A. Unfortunately, quantitative understanding of such high order correlation 
functions is far beyond our present knowledge. We will therefore need to make some further 
approximations in order to get a tractable expression for the pressure. 



1. Decoupling of the n-body correlations 

In order to simplify the problem of the n-body correlations, we invoke a so called 'linear' 
decoupling approximation, which attempts to describe the n-body correlation function in 
terms of n — 1 two body correlation functions: 

ff W (r U) r 23 , • ■ ■ , r n _ 1>n ) = g^\r l2 )g^ {v 2 ,) • . . g W (r n _ 1>ri ) (33) 

Here it should be understood that (^^(r^) stands for the pair correlation function of 
monomers of type A with monomers of type B in the multicomponent mixture of asso- 
ciating monomers, g ( - 2 \r 2 3) stands for the pair correlation function of monomers of type B 
with monomers of type C and so on. Still, these two body correlation functions may be 
quite difficult to obtain. In order to simplify the problem, consider one of these correlation 
functions, say, (^(r^), in the limit of zero density: 
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g {2 \v l2 - p = 0) = exp(-K(r 12 ) + $(n 2 )]/fcT) (34) 

It can be seen that, in this limit, the pair correlation function may be exactly expressed 
in terms of the pair correlation function of a reference system with no bonding potential, 

(2) 

9o ( r i2), times the Boltzmann factor of the bonding potential: 

^ (2) (r 12 ;p = 0) =^ 2) (r 12 ;p = 0)exp(-$(r 12 )/A;T) (35) 

Using the linear approximation to the n-body correlations and considering the above equa- 
tion to hold true at any density, the A integral is simplified considerably, giving: 

A = 5 71 - 1 (36) 

where S is defined as: 

6 = f gi 2) (r 12 ) exp(-$(r 12 )/£;T)d 3 r 12 (37) 

2. An equation of state for the n-mer in terms of the thermodynamics and structure of the 

monomers 

As a consequence of the two approximations given for the n-body correlations, we are 
now able to write down an expression for the pressure of the n-mer in terms of the properties 
of the reference fluid. Indeed, after setting a — 1, substitution of Eq. |3B| into Eq. |32] leads 
finally to the following result: 

^ = - [n - !][!+,_] (38) 

Note that as we are considering the limit of complete association, p is equal to the polymer 
number density. 

Adding the contribution of the reference system to the previous equation, we are now 
able to write down an equation for the compressibility factor, Z = p/kTp of the fluid of 
n-mer s: 
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nZ - (n — 1)[1 + p 



dhiS 



(39) 



n— mer 



dp 



where Z is the compressibility factor of the reference fluid, measured at the same monomer 
density as the n-mer fluid. This is a rather remarkable equation, as it gives the equation of 
state of a chain fluid from the properties of a fluid of monomers alone. Different versions 
of this equation will arise from the different theories available to describe the structure and 
thermodynamics of the fluid of monomers. In section IV we will consider two such theories 
in order to describe our model polymer. Let us recall at this point, however, that a simple, 
qualitative version of Eq.|39| may be obtained by simply considering that 5 does not depend 
on the density. In this way, the resulting equation does no longer depend on the structure 
of the reference fluid. Nikitin et a/.S have explored this equation using the van der Waals 
equation of state for Z and find the same qualitative behavior as is found in this work. 



The arguments we have put through in order to arrive at Eq. |3^ are rather physical and 
intuitive. On the other hand, Wertheim has developed a very general theory of association 
based on a re-summed cluster expansion, where the significance of each of the approximations 
is mathematically well understood. It is interesting to compare the results of this rather 
formal theory with the physically appealing description that we have used, largely based on 
the work of Zhou and StellS. 

In the extension of Wertheim's theory of association, the compressibility factor of the 
chain molecule is given as: 



F. Comparison with Wertheim's theory of association 




(40) 



where k is defined as: 




(41) 
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The limit of complete association requires that $ have an infinite well depth, so, within 
most of the range of the bonding potential, the Boltzmann factor is exceedingly bigger than 
1. Therefore, k and 5 become identical and Eq. ^0] is essentially equal to Eq. |39|. 

Before proceeding to the next section, let us first summarize the approximations invoked 
to obtain Eq. |SP|: 

1 . Assume that the free energy difference between the reference system and the completely 
associated system takes the form of the low density limit in all the density range. 

2. Decouple the n-body correlation function of the associating system into n—1 two body 
correlation functions through a linear approximation. 

3. Assume that the two body correlation function of the associating system is given in 
terms of the two body correlation function of a reference system with no bonding 
potential, as suggested by the exact low density limit. 

III. PREDICTIONS FOR THE SCALING LAWS OF THE CRITICAL 

PROPERTIES 

We start by assuming that the critical density does become small for large chain lengths, 
so that one can describe the equation of state in terms of a truncated virial expansion. 

l- = p + B 2 (T)p 2 + B 3 (T)p 3 (42) 

where p is the polymer number density. By applying the conditions for the critical point of 
pure fluids, i.e., 

we obtain a set of equations for the critical temperature and density@ 



B 2 (T C ) + JW 3 (T C ) = (44) 



3B 3 (T c ) Pc = 1 (45) 
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By making a Taylor expansion on powers of the density, the first and second virial 
coefficients predicted by Wertheim's equation are found to be: 



B 2 = n 2 (& 2 - ^a 2 ) 
B 3 = n 3 (b 3 - ^±a 3 



(46) 



where b 2 and 63 are the second and third virial coefficients of the reference fluid of non- 
bonded monomers, while a 2 and 03 are the zeroth and first order coefficients in a monomer 
density expansion of dlnd/dp. Of course, all these quantities are chain length independent. 

Now, in order to solve Eq. ^4] for the critical temperature we will need to linearize the 
virial coefficients with respect to the temperature. To do so, let us assume for the time 
being that there is a finite asymptotic critical temperature in the limit of infinite chain 
length, which we call 0, in analogy with the polymer + solvent case. We now make a series 
expansion of B 2 and B 3 in powers of AT = — T up to first order, and consider the limit 
of this expression for large n, leading to 

B 2 (T) = n 2 (C 2 - C' 2 AT) 
B 3 (T)=n 3 (C3-C3AT) 

where C 2 = b 2 (Q) — a 2 (0) and C 3 = 63(6) — 03(G) while C' 2 and C 3 are the corresponding 
derivatives with respect to temperature. Substitution of the linearized virial coefficients into 
the condition for the critical temperature leads to a quadratic equation for AT. Solving for 
this equation yields AT c (n), defined as 9 — T c (n): 

^ ~ I - ±w* i 12 ^ ~ 12C ^ + O " 2 7k ~ § k (48) 

This equation shows that AT c (n) must reach an asymptotic finite value, since the right hand 
side term should ultimately vanish for large n. The requirement for T c (n) to attain a finite 
asymptotic critical temperature equal to is then obeyed provided that C 2 vanishes. If we 
now notice that 

lim B 2 (&) =n 2 C 2 (49) 
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we arrive at the conclusion that indeed C 2 must vanish at the Boyle temperature of the 
infinitely long polymer, TgP, thus identifying the temperature with the Boyle temperature 
of the infinitely long polymer. From the definition of C 2 we see that this temperature is 
attained when the following condition is obeyed: 

b 2 (Tg) - a 2 (Tg) = (50) 

Note also that the leading terms of the expansion (Eq. ESI) are of order n _1//2 and just 
as predicted by the FH theory. 

The case of the critical polymer density is much simpler. Substitution of the expression 
for B 3 in the condition for the critical density shows that: 

p c (n) oc n- 3/2 (51) 

so that the critical mass density decreases with a power law proportional to n -1 / 2 , as pre- 
dicted by the FH theory. 

It is important to note that the above arguments apply regardless of the specific form 
in which the reference fluid (thermodynamics and structure) is described. In particular, the 
simplest implementation of TPT1, proposed by Nikitin et aZ@ considers 5 to be a constant. 
In such a case, a 2 is zero at all temperatures. However, Eq. [5(] shows that this simple version 
still predicts an asymptotic critical temperature which must obey the condition b%(T^) = 0. 
Obviously, this condition is obeyed for the Boyle temperature of the reference fluid. 

Another interesting issue is the apparent universality of the compressibility factor as 
predicted by the truncated virial expansion of Eq. [|2|. Indeed, substitution of this equation 
into the condition for the critical point shows that, apart from Eq. it must also hold 
that p c = —B 2 /(3B 3 ). Using this expression for the critical density in the linear term of Eq. 



42| and Eq. [|5| in the quadratic term, it is seen that both terms cancel each other exactly. 
Dividing the resulting expression for the pressure by the critical density (Eq. f|5]) then shows 
that: 

Z ^ = l^r=\ + W^ + - (52) 



Obviously, this result is independent on whatever assumption is made concerning the actual 
n dependence of the virial coefficients and shows that a finite asymptotic critical compress- 
ibility factor of about 1/3 is expected in the limit of infinite chain length, irrespective of the 
nature of the polymer. In the context of TPT1, a constant compressibility factor implies 
that the critical pressure must decrease as n~ 3 / 2 . 



IV. APPLICATION TO A POLYMER MODEL 

Let us consider a polymer model as that described in the previous section, with the 
reference fluid considered to be a truncated and shifted potential of the form: 



« (r) 



V LJ (r) - V LJ (r c ) r<r c 

(53) 

r > r r 



where r c = 2 • 2 1 / 6 and Vlj is the usual Lennard- Jones potential, 



a\ 12 



VWr)=4^ 7 j - W J (54) 

As to the bonding potential responsible for the connectivity between adjacent monomers, 
we will consider the FENE potential, defined in terms of Rq, the maximum displacement 
between monomers and k , a sort of elastic constant: 



$(r) 



-k Rl ln(l - 4, ) - E b < r < i? 

*° (55) 

r > Rq 



In order to ensure permanent connectivity of the n-mer, a constant, Ej,, which is (conceptu- 
ally) made infinitely large, is added to the actual FENE potential. 

In what follows, we will set k = 15e/a 2 and R = 1.5a and use the Lennard- Jones energy 
and range parameters as energy and length units, respectively. At liquid-like densities the 
most probable distance between non-bonded monomers is about 1.12a, which is bigger than 
the most probable distance 0.96a between bonded monomers. Note at this point that most 
of the previous applications of Wertheim's theory have been restricted to bonding potentials 
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of infinitely short range, allowing for a single possible bond length. To our knowledge, only 
once has the effect of a soft bonding potential been considered previously.^ 

In section II we related the equation of state of such a polymer fluid to the properties 
of the reference system of LJ monomers. What is now required is a theory for both the 
thermodynamics and the structure of the monomer fluid. We have obtained the required 
input from integral equation theory and thermodynamic perturbation theory. Let us consider 
each of them in turn. 



A. The RHNC integral equation theory 

In this approach, one attempts to calculate the exact pair correlation function, which is 
then used to evaluate the mechanical properties of the fluid. The Ornstein-Zernike equation 
relates the total pair correlation function, h(r) — g(r) — 1 to a short range direct correlation 
function c(r): 

h(r 12 ) = c(n 2 ) + p J h(r 13 )c(r 2 3)d 3 r 3 (56) 

Additionally, this integral equation must be provided with a closure that relates c(r) to 
h(r). We use the Reference Hyper-netted chain equation of Lado and Ashcroft.SH Although 
this set of equations can only be solved numerically and convergence is not a trivial matter, 
an efficient algorithm due to Labik and Malijevskyil makes the calculations affordable with 
a modest amount of CPU time. Once g(r) is known, the pressure of the fluid may be 
calculated using the standard relation^: 

V P 2 f d« 



k B T p -6^TJ r ^ 9{r) 4 " rdr (5?) 

Similarly, 5 may be calculated using Eq. [37|. Actually, solving the integral equation for each 
of the desired thermodynamic states may result rather cumbersome. In practice we solve 
the OZ+RHNC equation for several hundreds of state points and fit the pressure and S. 
Details of the procedure may be found in Appendix A. 
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B. The perturbation theory of Tang and Lu 



Perturbation theory was the first approach to give quantitative results for the thermo- 
dynamics of simple fluids at high density.ll'il Compared to integral equation theory, it gives 
similar results at high densities at a smaller computational cost, with the advantage that 
the free energy is obtained directly, without the need for thermodynamic integration. On 
the other hand, the traditional perturbation theories of Barker and Henderson@ and Weeks- 
Chandler- AndersenS are known to be rather poor at low densities because the underlying 
assumptions of these theories no longer hold true.0 Fortunately, Tang and Lu have pre- 
sented rather recently a second order perturbation theory for Lennard- Jones fluids which 
is very accurate both at low and high densitiesBH The success of this theory relies on a 
rather good description of the structure of the fluid, which is obtained from the OZ equation, 
supplemented by a simple closure known as the Mean Spherical Approximation. The use of 
this closure is very convenient because it has allowed to obtain a very good approximation 
to the actual free energy with a purely analytical equation. Furthermore, Tang and Lu 
have been able to obtain also analytic expressions for the pair correlation function using 
the MSA closure.00 With minor modifications we were able to extend this theory to our 
truncated and shifted Lennard- Jones potential and obtain a lengthy but analytic expression 
for the thermodynamics and structure of our reference fluid. Details of the implementation 
are explained in Appendix B. 

In what follows, we shall present the results of Wertheim's theory for our polymer model 
using both the RHNC and the MSA thermodynamic theories for the reference fluid. We will 
call each of the versions TPT1-RHNC and TPT1-MSA, following the original name for Eq. 
391 due to Wertheim.0 
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V. SIMULATION DETAILS 



In order to test the TPT1 theory, we have performed extensive computer simulations. 
Chain length n = 10 was chosen for a detailed comparison to the theory. We have calculated 
the pressure and chemical potential for 5 temperatures T = 1.68,2.5,3.0,4.0, and 5. The 
lowest value corresponds to a subcritical isotherm, while the highest value is above the G 
temperature. The pressure isotherms were evaluated from the virial in NVT Monte Carlo 
simulations The length L of the (cubic) simulation box was fixed to 18 a units. The density 
dependence of the chemical potential was calculated using grand-canonical Monte Carlo 
simulations. Chain conformations were sampled using local monomer displacements and 
slithering snake like motions. Particle insertions and deletions were performed following 
configurational bias grand-canonical acceptance rulesfl'Hil In order to obtain the equation 
of state we employ cycles of 25 local moves, 25 reptations and 10 CBGC moves. About 
40000 such cycles were performed so that at least a few thousand particle insertion-deletion 
attempts were accepted. The volume of the simulation box was chosen so that an average 
number of about 50 chain molecules was obtained. 

In order to compare the theory to the simulations, we must make sure that the chemical 
potentials are expressed with respect to the same reference state. In order to do so, we define 
the excess chemical potential to be the difference between /i and the chemical potential of 
an ideal gas of chains /z^ with the full intramolecular interactions but no intermolecular 
interactions. 

/i ex = /i - /i id with /i id = lnp — (n - 1) In C — In (Wo) (58) 

The first term denotes the translational entropy. Contributions due to the integration over 
the momenta are ignored, because they contribute equally to the reference system and the 
interacting polymer liquid. To determine the ideal gas contribution we construct chains 
according to the Rosenbluth procedure. The distance I between bonded neighbors is chosen 
according to its Boltzmann weight p(\l\) = 4:irl 2 exp(— (uq + <&)/kT)/C where C is the nor- 
malization constant. Wq denotes the Rosenbluth weight of the chains due to non-bonded 
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interactions, measured at zero density. Once the excess chemical potential has been ob- 
tained, it is compared to the excess chemical potential as predicted by the theory, which is 
evaluated using the standard thermodynamic relation: 

, . ex A ex 

^— = — + Z-1 (59) 

k B T Nk B T 1 ; 

The grand-canonical ensemble allows also for an accurate measurement of the phase 
diagram, because the order parameter (i.e., the density) is not conserved and density fluc- 
tuations are efficiently equilibrated. We monitor the probability distribution P(p) of the 
density. Close to two phase coexistence, the probability distribution is bimodal: one peak 
corresponds to the vapor, the other corresponds to the liquid. The coexistence chemical 
potential p C oex is fixed by the condition of equal weight in both peaks:^! 

/ dp P{p) = / dp P(p) with p* = dp pP{p) (60) 

JO Jp* Jo 

Far below the critical points the probability between the two peaks is very low, and we use 
a re-weighting scheme as to encourage the system to "tunnel" between the two phases. To 
this end we add a term kBTlnW(p) to the original Hamiltonian. Choosing W(p) ~ P(p) 
the system visits all densities with roughly equal probability. The probability distribution 
of the grand-canonical ensemble is obtained via re-weighting the distribution in the simula- 
tions Pmc according to P(p) = Pmc(p)W (p) . At very low temperatures the density of the 
liquid in coexistence with its vapor becomes very high and the configurational bias scheme 
becomes quite inefficient. Since the density of the vapor is very low, however, its pressure 
is vanishingly small. Hence, we employed NpT simulations at zero pressure to obtain the 
liquid density at coexistence. 

At the critical point the correlation length of density fluctuations diverges and universal 
behavior is expected. For finite chain length the unmixing transition exhibits 3D Ising 
universal behavior. We have located the critical point for chain length n — 1, 10, 20, 40, and 
60 by mapping the symmetrized order parameter distribution P sy m(p) = [P(p) +P(Pc~ p)]/2 
onto the universal scaling function of the 3D Ising model. This symmetrization reduces field 

23 



mixing corrections which are antisymmetric in p — p c to leading order. Normalizing P sym to 
unit variance and norm we eliminates all non-universal factors. The results of this mapping 
are presented in Fig.|2|, where we have used system sizes L = 11.3, 13.8, 18, 22.5, and 27 for 
chain length n — 1, 10,20,40, and 60, respectively. This method gives an accurate location 
of the critical temperature and density (finite size corrections to the critical density of the 
order L~^ 1 ~ a ^ v are neglected^). The locations of the critical points are collected in Tab. |. 

VI. RESULTS AND DISCUSSION 

Let us first examine the thermodynamic data for the chains of 10 monomers. Fig|J 
shows the predictions of TPT1 for several pressure isotherms (kT/e=5, 4, 3, 2.5 and 1.68) 
compared with simulation results. Both the RHNC and MSA versions of the theory are seen 
to give rather good estimates; at the highest temperatures, far above the estimated O point 
of our model (see below) as well as at the lowest, a subcritical isotherm. Overall, the RHNC 
version seems to describe the isotherms slightly better. 

Results for the excess chemical potential of the chains are shown in Fig.[|. The agreement 
is also quite satisfactory, though the results are slightly worse than for the pressure isotherms, 
specially at the lowest temperatures and densities. Indeed, the main assumption of the 
theory, that the local environment of a monomer in the polymer fluid is similar to that 
of the monomer fluid breaks down in the low density limit. The fluid is then made of 
isolated clusters of n monomers, rather than of single monomers uniformly distributed in 
space. Likewise, the theory is unable to describe the density dependence of the single chain 
internal energy and entropy. 

The liquid-vapor coexistence curve of the 10-mer as obtained from simulation and theory 
is shown in Fig.|5]. Both the RHNC and MSA versions overestimate the critical temperature 
as obtained from simulation by about 15%. Of course, this is expected for any classical 
theory. On the other hand, far away from the critical point, results from both versions of 
the theory are seen to yield fair agreement with simulation. The MSA version is somewhat 
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more convenient, however, because it allows to calculate the coexistence at low temperatures 
with no additional cost, while it becomes rather problematic to calculate the coexistence for 
the RHNC version below the reference fluid critical temperature. The reason for this is that 
the RHNC integral equation presents a region of no solutions below this point, so that the 
resulting equation of state is no longer defined inside the liquid-vapor envelope. 

We have also investigated the critical points of chains of 20, 40 and 60 n-mers, in an 
attempt to study the behavior of TPT1 for longer chains. Table | gives a summary of the 
simulation results, obtained by finite size scaling, together with predictions from TPT1- 
RHNC and TPT1-MSA. Both versions overestimate the critical temperatures by about 15% 
for all chain lengths studied. However, the MSA and RHNC predictions seem to converge 
as the chain length increases. On the other hand, the critical monomer densities are always 
underestimated, though the MSA version seems to give much better agreement than the 
RHNC version. In the latter theory the density decreases much too fast compared to the 
MC results. The overall behavior of the critical parameters is illustrated in Fig. [| where 
both T c and p c are plotted against n _1//2 , the predicted asymptotic scaling law for both of 
these properties. It is seen that for chain lengths up to 60 monomers, the critical properties 
are far from reaching their asymptotic behavior, so that the simulations do not allow as to 
asses unambiguously the predicted scaling laws. 

Although the calculation of the critical point of fluids larger than about 100 monomers 
by computer simulation becomes prohibitively expensive, we can estimate the O point of 
our polymer model by an analysis of the temperature dependence of the polymer extension. 
Fig. |7| shows a plot of the mean squared end to end distance divided by n — 1 as a function of 
temperature for various chain lengths. In the infinite chain length limit, the intercept of two 
such plots occurs at the O point of the polymer model. Extrapolation of the results gives as 
an estimate ~ 3.3. As to the theory, fitting the critical temperature predicted by TPT1- 
RHNC to a power law of the form T c = T c °° + bn' 1 ^ 2 + cn~ x in the range 10 2 to 10 7 gives 
T c °° = 3.44. On the other hand, by searching for the root in Eq. [50], we find that TPT1-MSA 
predicts T c °° = 3.14. Assuming that the O point is indeed the critical point of the infinitely 
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long chain, as suggested by the considerations of Section III, it would seem that TPT1 is 
capable of giving an excellent prediction for the point of the polymer, even though the 
actual prediction may vary somewhat depending on the theory used to describe the reference 
fluid. Remarkably, considering that 5 is density independent and using the simple van der 
Waals equation of state, Nikitin et al.0 have shown that TPT1 predicts an asymptotic 
critical temperature T c °° = yT c °. As T c °, the critical temperature of the reference fluid, 
is approximately 1.0 L J reduced units, we find that the simplest TPT1 approach already 
gives an excellent prediction for the G point of T c °° = 3.375. This is, however, somewhat 
fortuitous as it was shown in Section III that this TPTl-van-der- Waals approach of Nikitin 
actually predicts that is equal to the Boyle temperature of the monomer fluid, T^, which 
is about Tg = 2.58. Thus, the relatively good estimate of T c °° turns out to be a consequence 
of the over prediction of Tg implicit in the van der Waals equation of state. 

Recently, Chatterjee and SchweizeiS have analyzed the behavior of the critical point of 
infinite chain lengths using the PRISM theory. For two of the closures employed, the same 
behavior as that predicted by TPT1 is observed, at least concerning i) vanishing critical 
density and ii) finite critical temperature. The power laws are, however, different. The 
critical monomer density is predicted to vanish with a weaker dependence which may be 
either —1/3 or —1/4 depending on whether the RMPY/HTA or the MSA closures are used. 
The critical temperature is predicted to reach a finite critical value with the same exponents 
as the critical density. However, there was no a priori reason for choosing one closure over 
the others and several different trends could be obtained depending on the closure that was 
used. It is pleasing to see that TPT1 is able to give a unique conclusion, independent of the 
molecular theory used to describe the monomer fluid. 



VII. CONCLUSION 

In this paper we have used the formalism of Zhou and StellS to extend the TPT1 theory 
to polymers with variable bond-length. By using rigorous molecular theories for the reference 
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fluid of non-bonded monomers we have been able to explore two versions of the theory that 
allow to give a good description of the fluid without the need of any empirical data for the 
polymer. Comparison with numerical simulations show that both the RHNC version and 
the MSA version give good agreement with the simulation data. At low temperatures, the 
RHNC version seems to be more reliable, while at high temperatures there is apparently 
little difference. The MSA version is seen to give fair predictions for the critical points of 
the longer chains, with the advantage that it is almost analytic. The point of the model 
is predicted in very good agreement by the RHNC version, as well as by the MSA version. 

Concerning the critical behavior of long chains, it has been shown that TPT1 predicts 
an approach of the critical temperature to the point with a power law of n~ 1 ^ 2 . This is 
the correct behavior for the infinitely long polymer chain, as mean field behavior must be 
recovered in the infinitely long chain limit. The point has been shown to be the Boyle 
temperature of the infinitely long polymer while the critical mass density is predicted to 
vanish with a power law of n~ x l 2 . All of these predictions concerning the scaling behavior 
of the critical points of the pure monomer fluid are seen to agree exactly with the mean 
field (Flory-Huggins) predictions for the critical behavior of polymer+solvent mixtures. 
This gives further support to the idea that pure polymer equations of state may be used 
as effective equations of state for the polymer+continuum-solvent system and vice versa, 
polymer+solvent equations of state may be used for the polymer+vacuum case; a formal 
prove of this intuitively appealing idea is, however, difficult. 

In a subsequent papei@ we use the implementation of TPT1 proposed here to describe 
the equation of state, together with a self consistent field theory to study the surface and 
interfacial properties of a polymer+solvent system. 
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VIII. APPENDIX A: FITTING THE RHNC DATA FOR THE MONOMER FLUID 

Rather than solving the integral equation of each of the state points desired, which 
would be rather expensive, we followed a similar approach as that used by Johnson et 
a/.il to describe the Lennard- Jones fluid. We solved the RHNC integral equation for a 
set of 756 states in the range 1.1 < T < 6.0 and < p < 0.85, we calculated the resulting 
compressibility factor (see Eq. [57]) and then fit the data to a Modified-Benedict- Webb- Rubin 
equation of state, given by: 

pkT(Z - 1) = E a iP i+1 + FJ2 hp 2i+1 (61) 

i=l i=l 

where F = exp(— 7p 2 ), while the exact form of the and hi coefficients is given in Table 0. 
Once the fit to Z is performed, the free energy may be determined from: 

A/N = lt — + £&A (62) 

i=l 1 i=l 

The Gi coefficients obey the following recursive relation: 

Gi = _ *>*'-> -ff- DC- (63) 

with the first term given by G\ = (1 — F)/(2j). The parameters obtained for the fit are 
collected in Table |TJ. 

Contrary to the approach of Johnson et at., we calculate the thermodynamics using the 
RHNC theory, rather than computer simulations. In this way we are able to save several 
orders of magnitude of CPU time. The effective diameter required in the RHNC equation 
is determined as suggested by Lado et al.^ while the effective hard sphere bridge function 
is calculated from the parameterization of Labik and Malijevskyffl 
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Also required is the associating strength 5 of the fluid. This is obtained by solving 
equation |37] for each of the 756 state points. Rather than fitting 5, which is a difficult task, 
as it varies several orders of magnitude in the range 1.1 < T < 6.0, we fit the dimensionless 
ratio 8(p)/5(p = 0). The actual functional form used is: 

5 5 

5/5, = 1 + E E ^P l T {1 - ]) (64) 

i=l j=l 

The parameters for this fit are gathered in Table |TV| 

IX. APPENDIX B: THE PERTURBATION THEORY OF TANG AND LU 

In order to describe the thermodynamics of the reference fluid we perform a Barker- 
Henderson decomposition of the monomer fluid pair potential such that uq is described in 
terms of a repulsive reference potential u> re f, which is made of the positive region of uq and 



a perturbation, w per , which is made of the negative part of the potentialHE. 

Uo(r) r < to r < to 

w re f(r) = \ and w pcr (r) = u (r) - w ie{ (r) = I (65) 

r > to \ uq{t) r > to 

where t = 1.0013 defines the value of r where Uq becomes negative (recall that we are 
considering a cut and shifted potential). We now couple the perturbation potential to the 
reference potential with a coupling parameter, A, so that the actual potential is recovered 
for A = 1: 

w(r; A) = w rci (r) + Xw per (r) (66) 

At this stage we recall the fundamental functional expression that relates the Helmholtz free 
energy with the radial distribution function,il 

5 A 1 

-Npg(r;X) (67) 



5w(r;X) 2 

Integration of this equation following the rules of functional calculus,il leads to an expression 
relating the free energy of the monomer fluid with that of the reference fluid: 
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A _ Aref = 1 I I ^ A)w pcr 4vrr 2 drdA (68) 

iV iV Z JX=0 J to 

It is then assumed that the radial distribution function may be expanded as a series in 
powers of A of the form g(r; A) = go(r) + gi(r)X + .... Truncation of the series to first order 
then yields: 

A A T f 1 f°° 1 f°° 

— -jp = -p g (r)w peT (r)4iir 2 dr + -p gi{r)w pcv {r)A^r 2 dr (69) 

Following Barker and Hendersonil we now choose to describe the reference potential by 
means of a hard sphere fluid of appropriate hard sphere diameter, d. This choice is justified 
because the reference potential is made essentially of the repulsive part of the potential. In 
this way, j3A re {/N and go may be considered to be the free energy and radial distribution 
function of an effective hard sphere fluid, while g\ may be solved using the MSA closured 
However, the integrals appearing in the previous equation are still quite tedious to calculate 
and a further approximation allows to get fully analytical results. This is done by fitting the 
actual monomer potential, u$ to a Two Yukawa potential, following the procedure of ref.il: 

u (r) TY = -he— r + he 6 — (70) 

Actually, the fit needs to be performed only for values of r greater than d and the resulting 
function (with k Q = 2.4405, Z\ = 3.492456 and z<± = 13.109857) is virtually identical to the 
true potential. However, substitution of Eq. ^ in the first and second order contributions 
of Eq. |69| and solving g% for u TY rather than for Uo, a very accurate expression for Eq. ^ 
may be obtained which is fully analytic. 

This is done by rearranging Eq. |69| into the form: 

— = a + a 1 + a 2 (71) 

where is given by the Carnahan-Starling equation of state: 

ao = JT^W (72) 

while 
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I poo I poo pd 

fl i = nP / (9o - l)wJj47rr 2 dr + -p w per 47rr 2 dr + g (d) / w per 47rr 2 dr (73) 

Z Jd Z Jd Jta 



and 



0-2 = ~P I gi(r)vjpe r 4irr 2 dr + -pgi(d) / w per 4irr 2 dr 

4 Jd 4 Jta 



(74) 



In the last two equations it is understood that wT^ is the perturbation potential when 
expressed as in Eq. [TCj; the actual w peT potential is used where ever possible; g — 1 and 
gi = beyond the cutoff distance of the potential while go and gi are considered to remain 
constant in the range [d, ta] . The only difference between the above expressions and those 
obtained by Tang et a/iiH for the true Lennard- Jones potential are found in the trivial 
integrals of the form / u>i47rr 2 dr because both the integration limits and the perturbation 
potential differ. Solving a± and a 2 yields: 



a 1 



I2rjpe 

~1F~ 



L(z\d) 



1 + Zid s 
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-h 



L(z 2 d) 
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(75) 
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where r] = %/6d 3 p and ki = fc e CT ■ The Barker- Henderson diameter, defined as 



(76) 



d= (1 - e-^W 
J 



)dr 



(77) 



is 



parameterized using the formula proposed in RefS 
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d/a = 2 1 / 6 



1+1 + 



T + c 2 T 2 + c 3 T' 

Cl 



4\ 1/2' 



-1/6 



where d = 1.150167, c 2 = -0.046498, and c 3 = 0.0004477054. 
On the other hand, Q is defined as: 

S(t) + UriL^e-t 



where 



Q(t) = 



(78) 



(79) 



S(t) = (1 - r/) 2 r + 677(1 - r/)t 2 + 18r]H - 12r/(l + 2?y) 



(80) 



and 



L{t) = (1 + rj/2)t + 1 + 2r) 



(81) 



The reference radial distribution at contact is given by: 



(82) 



while W cs is to constant factors, the definite integral of the perturbation potential, w per: : 



Wc(x,y) 



a 



x 



a 



x 



x\ 3 ( 3 



a 



(83) 



The compressibility factor may be determined by density differentiation of the free energy, 
which leads to: 



pV 
NkT 



Zq + Z\ + Z2 



(84) 



where 



1 + rj + rf — rf 



(85) 
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Z x = at 
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(86) 
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6(1 - rj)t 2 + 36r/t - 12(1 + 577) + 12[(1 + 2r;)t + 1 + 5r7]e" 



(87) 



(88) 



(89) 



(1 -r])H 3 

In order to calculate the associating strength, 6, the radial distribution function could 
have been assumed to be g = go + g±, which is already a rather good approximation.il How- 
ever, we use the SEXP (simplified Exponential) approximation which considerably improves 
the estimate of g around a with no additional information.il According to this approxima- 
tion, g = goe 91 . Once g is known, 5 is calculated directly by invoking Eq. |S7|. The integral 
must be performed numerically but would have been analytic if the bond length was held 
fixed. 
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FIGURES 

FIG. 1. Relating the structure of the fluid to the excess chemical potential of the components 
by means of a thermodynamic cycle. 

FIG. 2. Mapping of the probability distribution of the density onto the universal distribution 
of the 3D Ising universality class (line). The chain length n and the estimate for the critical 
temperature are indicated in the key. 

FIG. 3. Pressure against monomer density for chains of 10 monomers. All quantities are given 
in Lennard-Jones reduced units. Symbols are NVT simulation data while lines are predictions from 
TPT1; full line, RHNC version; dashed line, MSA version. From top to bottom, pressure isotherms 
at T = 5, 4, 3, 2.5 and 1.68 reduced LJ units. 

FIG. 4. Excess chemical potential against monomer density for chains of 10 monomers. All 
quantities are given in Lennard-Jones reduced units. Symbols are Grand Canonical simulation 
data, while lines are predictions from TPT1; full line, RHNC version; dashed line, MSA version. 
From top to bottom, chemical potential isotherms at T = 5,4,3,2.5 and 1.68 LJ reduced units. 

FIG. 5. Liquid-vapor coexistence curves of a 10-mer as obtained from grandcanonical simula- 
tions (solid lines) and NpT simulations at p = (diamonds), compared with TPT1-RHNC (dashed 
line) and TPT1-MSA (dotted line). The filled circle presents the critical point as extracted from 
finite size scaling of the MC data. The open circle and the open square denote the critical point 
of the TPT1-RHNC and the TPT1-MSA, respectively. The simulation results and TPT1-MSA 
calculations for monomers (n = 1) are also included. 

FIG. 6. Critical temperatures (a) and critical monomer densities (b) in the MC simulations 
and the perturbation theory. In panel (a) the temperatures of the intersections of R^{T)/(n — 1) 
for neighboring chain lengths are also included. For n — > oo the values tend to the temperature. 
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FIG. 7. Temperature dependence of end-to-end distance of a single chain in the vicinity of 
the temperature. The crossings of the ratios R^{T)/{n — 1) for neighboring chain lengths are 
indicated by arrows. The crossing points converge to the G temperature in the limit of infinite 
chain length n — ► oo. 
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TABLES 



n 


T c (MC) 


T c (RHNC) 


T c (MSA) 


1 


1.00 


1.02 


1.11 


10 


1.98 


2.27 


2.36 


20 


2.214 


2.56 


2.62 


40 


2.396 


2.79 


2.81 


60 


2.485 


2.90 


2.88 


oo 


~ 3.3 


3.44 


3.14 


n 


np c (MC) 


n Pc (RHNC) 


n Pc (MSA) 


1 


0.321 


0.376 


0.323 


10 


0.245 


0.207 


0.217 


20 


0.206 


0.145 


0.184 


40 


0.172 


0.108 


0.150 


60 


0.1523 


0.091 


0.140 




TABLE I. Critical temperature, T c and critical 


monomer density, np c 


as obtained from simu- 



lation (MC) and from TPT1 with either the RHNC version or the MSA version for the structure 
and thermodynamics of the reference fluid 
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1 x x T + x 2v /T + x 3 + Xi/T + x 5 /T 2 x 20 /T 2 + X21/T 3 


2 x 6 T + x 7 + x 8 /T + x 9 /T 2 


X22/T 2 + x 23 /T 4 


3 x w T + x u + xn/T 


X24/T 2 + x 25 /T 3 


4 x 13 


x 26 /T 2 + x 27 /T* 


5 x u /T + x 15 /T 2 


X28/T 2 + x 29 /T 3 


6 xi 6 /T 


x 30 /T 2 + x 31 /T 3 + x 32 /T 4 


7 x 17 /T + x 18 /T 2 




8 x 19 /T 2 





TABLE II. The ctj and hi temperature dependent coefficients of the BWR equation of state 



(61). The Xj are adjustable parameters whose actual value are given in Table 111 
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1 


Xj 


i 


X i 


1 


.421192000D+00 


17 


-.107917493D+03 


2 


.476645300D+01 


18 


.946452266D+03 


3 


-.915420500D+01 


19 


-.661602970D+03 


4 


.161011000D+01 


20 


.117245702D+03 


5 


-.152328700D+01 


21 


-.348234700D+01 


6 


.205245400D+01 


22 


.415445097D+03 


7 


-.348800300D+01 


23 


.325092260D+02 


8 


.459058600D+01 


24 


.136021606D+04 


9 


-.111109576D+03 


25 


-.672466920D+03 


10 


-.683045000D+00 


26 


-.114228784D+03 


11 


.106678670D+02 


27 


.602317138D+03 


12 


-.243251130D+02 


28 


.185990661D+03 


13 


.193579400D+02 


29 


-.116953259D+04 


14 


-.196790929D+03 


30 


.607836000D+00 


15 


-.172144369D+03 


31 


.656689000D+00 


16 


.262199061D+03 


32 


.588415000D+00 



TABLE III. Parameters for the fit of the RHNC data to the Modified Benedict- Webb-Rubin 



equation of state. Notation as in paperE3. The nonlinear parameter is set to 7 = 3 
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i 


1=1 


i=2 


i=3 


1=4 

J 


i=5 


1 


.496554000D+00 


.129871000D+01 


-.595011800D+01 


.860280100D+01 


-.452759900D+01 


2 


.101561200D+01 


-.967002900D+01 


.677884460D+02 


-.950135670D+02 


.471191210D+02 


3 


-.604253400D+01 


.411972640D+02 


-.252396610D+03 


.376561440D+03 


-.178249974D+03 


4 


.622698400D+01 


-.680522730D+02 


.391425109D+03 


-.605492785D+03 


.292582291D+03 


5 


-.691995000D+00 


.357118100D+02 


-.218469635D+03 


.348511635D+03 


-.171957868D+03 



TABLE IV. Parameters a%j for the fit to S/Sq (Eq. [34|) as obtained from the RHNC integral 
equation 
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